Watersheds
2.0 Overview
The following compares two R-based workflows for watershed delineation and hydrological network extraction from digital elevation models (DEMs). We implement parallel processing pipelines using flowdem (Martinsen, 2023) and whitebox (Wu & Brown, 2025) packages to evaluate algorithmic differences in depression handling, flow routing, and stream network generation. Both workflows leverage high-performance geospatial libraries including RichDEM (Barnes, 2016), WhiteboxTools, GDAL 3.11 (contributors, 2020), and the r-spatial ecosystem (Hijmans, 2025; Pebesma & Bivand, 2023a, 2023b).
The comparative approach serves two purposes: (1) validating watershed boundaries through algorithmic consensus, and (2) identifying optimal processing strategies for challenging terrain features such as endorheic basins and artificial depressions. Performance benchmarks and visual comparisons guide selection between rapid exploratory analysis (flowdem) and fine-tuned parametric control (whitebox).
Environment Setup
pacman::p_load(
"bslib",
"cli", "cols4all", "covr", "cowplot",
"dendextend", "digest", "DiagrammeR",
"dtwclust", "downlit",
"exactextractr", "elevatr",
"FNN", "future", "flowdem",
"gdalUtilities", "geojsonsf", "geos", "geodata",
"ggplot2", "ggstats","ggspatial", "ggmap",
"ggplotify", "ggpubr", "ggrepel", "giscoR",
"hdf5r", "httr", "httr2", "htmltools",
"jsonlite",
"leafem", "leaflet.providers", "libgeos",
"luz", "lwgeom", "leaflet", "leafgl",
"mapedit", "mapview", "maptiles",
"methods", "mgcv", "MPI",
"ncdf4", "nnet",
"openxlsx",
"parallel", "plotly", "proj4", "PROJ", "progress", "purrr",
"randomForest", "rasterVis", "raster",
"rayshader", "rayvertex",
"RColorBrewer", "rgl", "rmapshaper", "rsconnect",
"RStoolbox", "rts", "rgrass",
"s2", "sf", "scales", "spdep", "stars",
"stringr", "supercells",
"terra", "terrainr", "testthat", "traudem", "taudem",
"tidyverse", "tidyterra", "tools",
"tmap", "tmaptools", "terrainr",
"whitebox", "xgboost"
)2.1 Declare AOIs
crs_master = sf::st_crs("EPSG:3857")
country = giscoR::gisco_get_countries(
country = "Malawi", resolution = "3") |>
sf::st_cast() |> sf::st_transform(crs_master)
lake = sf::st_read("../assets/inputs/lakes_site.shp") |>
sf::st_cast() |> sf::st_transform(crs_master)
NA Reading layer `lakes_site' from data source
NA `/Users/seamus/repos/map-templates/assets/inputs/lakes_site.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 1 feature and 22 fields
NA Geometry type: POLYGON
NA Dimension: XY
NA Bounding box: xmin: 35.6 ymin: -15.5 xmax: 35.9 ymax: -15
NA Geodetic CRS: WGS 84
bbox = lake |>
sf::st_buffer(dist = 60000) |>
sf::st_bbox() |>
sf::st_as_sfc() |>
sf::st_sf()
# Interactive map mode: "view"
tmap::tmap_mode("view")
tmap::tm_shape(country) + tmap::tm_borders(lwd=1, col= "green") +
tmap::tm_shape(bbox) + tmap::tm_borders(lwd=2, col= "orange") +
tmap::tm_shape(lake) + tmap::tm_borders(lwd=2, col= "blue") +
tmap::tm_basemap("Esri.WorldImagery")Figure 2a: Interactive map illustrating layout of area of interest polygons (AOI)
2.2 Download DEM
We acquired elevation data using the elevatr package (Hollister et al., 2023) that accesses via Amazon Web Services Terrain Tiles and the Open Topography Global DEM API collections of global digital elevation models including SRTMGL3, SRTMGL1, AW3D30, and SRTM15Plus. Depending on source, OSM’s Slippy Tiling syntax1 is used where zoom represents number of vertical tiles counting from top-left corner in a Pseudo-Mercator grid (EPSG:3857).
At our study latitude, zoom=11 returns approximately 3-arc-second resolution (~76m at the equator). We then resample to a standardized 100m grid to facilitate downstream metrics computed at landscape scale. The crs_master variable ensures coordinate reference system consistency throughout subsequent operations.
# z = 12: 1-Arc Second Resolution
# z = 11: 3-Arc Second Resolution
# z = 10: 5-Arc Second Resolution
dem = elevatr::get_elev_raster(bbox, z=10, clip="locations") |> terra::rast()
names(dem) = "elevation"
dem_100m = stars::st_warp(
stars::st_as_stars(dem),
cellsize=100, crs=sf::st_crs(crs_master)) |>
terra::rast()
# Outputs
terra::writeRaster(dem_100m, overwrite = TRUE,
"../assets/TIF/dem_chilwa_00_raw.tif")2.3 Hydraulic Conditioning
Both flowdem and whitebox workflows follow the standard hydrological processing sequence: (1) hydraulic conditioning through depression handling, (2) flow direction and accumulation calculation, and (4) stream network extraction. However, they differ fundamentally in algorithmic implementation and parametric control.
flowdem integrates RichDEM operations using priority-flood algorithms (Barnes et al., 2014) with minimal parameters that are optimized for computational efficiency and exploratory analysis. Meanwhile, whitebox provides granular control through extensive tools, enabling parameter tuning for complex terrain such as karst formations or anthropogenic modifications (Hasan et al., 2021).
2.3.1 Depression Breaching
DEMs contain spurious depressions from acquisition artifacts, interpolation errors, and measurement limitations (Lindsay, 2016). These disrupt flow connectivity essential for watershed delineation. Hydraulic conditioning removes artifacts while preserving genuine topographic features like Lake Chilwa’s endorheic basin.
The flowdem workflow implements Barnes et al.’s priority-flood algorithm (Barnes et al., 2014) via the RichDEM library, which processes depressions via optimized queue-based operations. The breach() function carves minimal-depth channels through topographic barriers, improving efficiency with in-memory operations.
The whitebox workflow provides additional parameter controls with which we add a breaching constraint (max_depth=3m) and a slope component (flat_increment=0.001) to ensure monotonic descent along breach channels. This parameter-controlled approach enables selective breaching based on depression characteristics. The flat_increment parameter performs similar function as epsilon used in subsequent flowdem::breach operations, only with greater manual controls of the slope gradient (Lindsay, 2016).
Performance benchmarking reveals whitebox breaching completes in 2.3 seconds versus 8.8 seconds for flowdem, likely reflecting differences between compiled C++ (WhiteboxTools) versus R-wrapped operations (RichDEM). However, flowdem’s approach requires fewer parameters, facilitating rapid exploratory workflows.
# Compare breaching algorithms
time_wb <- system.time({
whitebox::wbt_breach_depressions(
dem = "dem_chilwa_00_raw.tif",
output = "dem_chilwa_11_breached_flat.tif",
wd = "../assets/TIF/",
max_depth = 3,
flat_increment = 0.001)
})
time_fd <- system.time({
dem_breach_fd <- flowdem::breach(
terra::rast("../assets/TIF/dem_chilwa_00_raw.tif"))
terra::writeRaster(dem_breach_fd, overwrite = TRUE,
"../assets/TIF/dem_chilwa_01_breached.tif")
})
# Compare runtimes
cat("\n=== Runtime Comparison ===\n")
NA
NA === Runtime Comparison ===
cat("Whitebox breach:\n")
NA Whitebox breach:
cat(sprintf(" Elapsed time: %.3f seconds\n", time_wb["elapsed"]))
NA Elapsed time: 2.344 seconds
cat(sprintf(" User time: %.3f seconds\n", time_wb["user.self"]))
NA User time: 0.014 seconds
cat(sprintf(" System time: %.3f seconds\n\n", time_wb["sys.self"]))
NA System time: 0.009 seconds
cat("flowdem breach:\n")
NA flowdem breach:
cat(sprintf(" Elapsed time: %.3f seconds\n", time_fd["elapsed"]))
NA Elapsed time: 8.772 seconds
cat(sprintf(" User time: %.3f seconds\n", time_fd["user.self"]))
NA User time: 8.577 seconds
cat(sprintf(" System time: %.3f seconds\n\n", time_fd["sys.self"]))
NA System time: 0.181 seconds
# Compare spatially
dem_100m = terra::rast("../assets/TIF/dem_chilwa_00_raw.tif")
dem_breach_fd = terra::rast("../assets/TIF/dem_chilwa_01_breached.tif")
dem_breach_wb = terra::rast("../assets/TIF/dem_chilwa_11_breached_flat.tif")
dem_breach_diff_fd <- dem_100m - dem_breach_fd
dem_breach_diff_wb <- dem_100m - dem_breach_wb
dem_breach_diff_fd[dem_breach_diff_fd == 0] <- NA
dem_breach_diff_wb[dem_breach_diff_wb == 0] <- NA
tmap::tmap_mode("view")
tmap::tm_shape(dem_breach_diff_fd) + tmap::tm_raster(values="Blues", title = "Change (m)",
breaks = c(-400, -100, -10, 0, 10, 100, 400),
labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
# tmap::tm_compass(color.dark="gray60,position=c("RIGHT", "top")) +
# tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_title("Breaching Effects (Barnes, 2016)", size=.8) -> tm_breach_fd
tmap::tm_shape(dem_breach_diff_wb) + tmap::tm_raster(values="Blues", title = "Change (m)",
breaks = c(-400, -100, -10, 0, 10, 100, 400),
labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
# tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
# tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_title("Breaching Effects (Lindsay, 2016a)", size=.8) -> tm_breach_wb
tmap::tmap_arrange(tm_breach_fd, tm_breach_wb, nrow=1)
Figure 2b: Maps comparing effects of breaching algorithms (Barnes, 2016; Lindsay, 2016). Both algorithms target similar depression features but with different modification magnitudes controlled by respective parameters.
Table 2a: Runtime comparison of breaching algorithms applied to Lake Chilwa watershed DEM (100m resolution, ~1.3M cells)
| Metric | Whitebox | flowdem | Performance Ratio |
|---|---|---|---|
| Elapsed time | 2.30s | 8.83s | 3.8× faster |
| User CPU time | 0.01s | 8.59s | 859× faster |
| System time | 0.01s | 0.21s | 21× faster |
The performance advantage of whitebox likely reflects architectural differences: WhiteboxTools implements breaching in compiled C++ with optimized spatial indexing, while flowdem wraps RichDEM operations through R interfaces that incur marshalling overhead. The high user CPU time in flowdem (97% of elapsed time) suggests computational delays in the R-wrapped priority-flood implementation, whereas whitebox delegates processing to external compiled binaries with minimal R overhead (0.4% of elapsed time).
Despite longer execution times, flowdem’s parameter-free interface eliminates decisions about depth thresholds and flat increments, potentially reducing workflow development time for exploratory analyses. For production workflows requiring repeated processing or large spatial extents (>10M cells), whitebox’s computational efficiency becomes increasingly critical.
Both algorithms targeted similar depression features across the study area (Figure 2b), though with different modification patterns. Whitebox’s max_depth constraint preserved deeper basin structures (>3m) while breaching shallow artifacts, producing elevation changes ranging from -10m to +10m in modified areas. Flowdem’s uniform breaching approach modified a broader spatial extent without depth discrimination, resulting in elevation changes up to ±100m in some locations. The hybrid nature of Lindsay’s algorithm, which combines selective breaching with subsequent filling, highlights whitebox’s more conservative modifications. In contrast, the priority-flood approach in RichDEM implements aggressive breaching to ensure complete depression removal (Barnes et al., 2014), which may be undesirable when preserving hydrological features like endorheic basins.
2.3.2 Depression Filling
Following breaching, remaining flat areas require treatment to ensure continuous flow paths. This step is particularly critical for endorheic basins like Lake Chilwa, where water accumulates through inflow and evaporation without surface outflow. Standard watershed algorithms typically assume exorheic drainage systems with defined outlet points (O’Callaghan & Mark, 1984; Wu et al., 2015). For closed basins, depression filling must preserve the terminal sink while ensuring flow connectivity throughout the contributing watershed.
flowdem implementation uses epsilon-gradient filling via fill(epsilon=T). This function adds infinitesimal elevation increments (typically 10⁻⁶ to 10⁻⁴ m) to flat surfaces, ensuring strict monotonic descent along flow paths (Barnes et al., 2014). The epsilon value also prevents ambiguous flow directions in perfectly flat terrain common to marshland environments surrounding endorheic lakes. Critically, the algorithm preserves major depressions like Lake Chilwa while establishing flow pathways toward this terminal sink.
whitebox implementation applies the scan-line algorithm (Wang & Liu, 2006) via wbt_fill_depressions_wang_and_liu() function. This method uses an efficient line-by-line approach optimized for low-relief terrain, filling remaining depressions after the initial breaching operation. When combined with conservative breaching (max_depth=3m), this two-step approach maintains the integrity of genuine hydrological features while removing artifacts that would fragment the drainage network.
The challenge for endorheic systems is distinguishing between genuine terminal sinks requiring preservation and spurious depressions requiring removal. Our approach addresses this through:
- Conservative breaching: Using depth constraints to preserve deep basins while removing shallow artifacts
- Selective filling: Applying fill operations to remaining flat areas without eliminating the primary drainage terminus
- Flow direction validation: Ensuring all cells establish flow paths toward the identified terminal sink
This conditioning strategy enables subsequent flow accumulation to correctly identify Lake Chilwa as the watershed’s convergence point rather than fragmenting drainage into multiple isolated depressions.
# Compare filling algorithms
time_wb <- system.time({
whitebox::wbt_fill_depressions_wang_and_liu(
dem = "dem_chilwa_11_breached_flat.tif",
output = "dem_chilwa_12_filled_wang.tif",
wd = "../assets/TIF/")
})
time_fd <- system.time({
dem_fill_fd <- flowdem::fill(
terra::rast("../assets/TIF/dem_chilwa_00_raw.tif"), epsilon=T)
terra::writeRaster(dem_fill_fd, overwrite = T,
"../assets/TIF/dem_chilwa_02_filled.tif")
})
# Compare runtimes
cat("\n=== Runtime Comparison ===\n")
NA
NA === Runtime Comparison ===
cat("Whitebox fill:\n")
NA Whitebox fill:
cat(sprintf(" Elapsed time: %.3f seconds\n", time_wb["elapsed"]))
NA Elapsed time: 1.930 seconds
cat(sprintf(" User time: %.3f seconds\n", time_wb["user.self"]))
NA User time: 0.008 seconds
cat(sprintf(" System time: %.3f seconds\n\n", time_wb["sys.self"]))
NA System time: 0.006 seconds
cat("flowdem fill:\n")
NA flowdem fill:
cat(sprintf(" Elapsed time: %.3f seconds\n", time_fd["elapsed"]))
NA Elapsed time: 3.411 seconds
cat(sprintf(" User time: %.3f seconds\n", time_fd["user.self"]))
NA User time: 3.362 seconds
cat(sprintf(" System time: %.3f seconds\n\n", time_fd["sys.self"]))
NA System time: 0.044 seconds
# Compare spatially
dem_100m = terra::rast("../assets/TIF/dem_chilwa_00_raw.tif")
dem_fill_fd = terra::rast("../assets/TIF/dem_chilwa_02_filled.tif")
dem_fill_wb = terra::rast("../assets/TIF/dem_chilwa_12_filled_wang.tif")
dem_fill_diff_fd <- dem_100m - dem_fill_fd
dem_fill_diff_wb <- dem_100m - dem_fill_wb
dem_fill_diff_fd[dem_fill_diff_fd == 0] <- NA
dem_fill_diff_wb[dem_fill_diff_wb == 0] <- NA
tmap::tmap_mode("view")
tmap::tm_shape(dem_fill_diff_fd) + tmap::tm_raster(
palette = "pi_yg",, title = "Change (m)", midpoint = 0,
breaks = c(-400, -100, -10, 0, 10, 100, 400),
labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
# tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
# tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_title("Filling Effects (Barnes, 2016)", size=.8) -> tm_fill_fd
tmap::tm_shape(dem_fill_diff_wb) + tmap::tm_raster(
palette = "pi_yg",, title = "Change (m)", midpoint = 0,
breaks = c(-400, -100, -10, 0, 10, 100, 400),
labels = c("-400 to -100", "-100 to -10", "-10 to 0", "0 to 10", "10 to 100", ">100")) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
# tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
# tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_title("Filling Effects (Wang & Liu, 2006)", size=.8) -> tm_fill_wb
tmap::tmap_arrange(tm_fill_fd, tm_fill_wb, nrow=1)
Figure 2c: Interactive maps comparing effects of filling algorithms (Lindsay, 2016; Wang & Liu, 2006).
Depression filling showed similar performance patterns to breaching (Table 2a). Whitebox completed filling in 2.0 seconds compared to 3.4 seconds for flowdem, representing a 1.7× speedup. Both algorithms maintained low system overhead, though flowdem again showed higher user CPU consumption.
Table 2b: Runtime comparison of filling algorithms
| Metric | Whitebox | flowdem | Performance Ratio |
|---|---|---|---|
| Elapsed time | 2.01s | 3.42s | 1.7× faster |
| User CPU time | 0.01s | 3.36s | 336× faster |
| System time | 0.01s | 0.05s | 5× faster |
The more modest performance difference compared to breaching (1.7× vs 3.8×) suggests Wang and Liu’s scan-line algorithm and RichDEM’s epsilon-gradient method have comparable computational complexity, with efficiency differences primarily attributable to implementation architecture rather than algorithmic approaches.
Both filling algorithms effectively smoothed flat areas while preserving Lake Chilwa’s endorheic basin structure (Figure 2c). Whitebox modifications concentrated along identified flat areas with conservative elevation adjustments (-10m to +10m), while flowdem produced more extensive modifications (up to ±100m) across broader spatial extents. The resulting elevation models from both methods achieved hydrological connectivity, though through different numerical pathways.
2.4 Flow Accumulation
Flow routing was calculated using the D8 method applied in the flowdem::dirs() and whitebox::wbt_d8_pointer() functions. The D8 single-flow-direction algorithm (O’Callaghan & Mark, 1984) routes flow from each cell to its steepest downslope neighbor among eight adjacent cells. Flow direction pointers are then used to compute flow accumulation values, which represent the number of upstream cells contributing flow to each location. At 100m resolution, each cell represents 0.01 km² (1 hectare) of contributing area, enabling direct conversion between cell counts and drainage area in km².
Flow accumulation calculations reveal critical differences between workflows. Whitebox achieves maximum accumulation of 1,152,591 cells (11,526 km²), indicating successful watershed-scale convergence to Lake Chilwa. In contrast, flowdem produces maximum accumulation of only 5,017 cells (50.2 km²), suggesting potential flow routing issues or incomplete drainage network development. This 230× difference warrants investigation, as both workflows used identical D8 algorithms and spatial resolution.
Runtime performance shows comparable efficiency: whitebox completes flow direction and accumulation in 2.0 seconds versus 2.8 seconds for flowdem, indicating both implementations efficiently handle D8 calculations at 100m resolution across the study extent.
time_wb <- system.time({
whitebox::wbt_d8_pointer(
dem = "dem_chilwa_12_filled_wang.tif",
output = "dem_chilwa_13_flow_direction_D8.tif",
wd = "../assets/TIF/")
whitebox::wbt_d8_flow_accumulation(
input = "dem_chilwa_13_flow_direction_D8.tif",
output= "dem_chilwa_14_flow_accumulation_D8.tif",
wd = "../assets/TIF/", pntr=T)
})
time_fd <- system.time({
dem_dir_fd <- flowdem::dirs(
terra::rast("../assets/TIF/dem_chilwa_02_filled.tif"), mode="d8")
dem_acc_fd <- flowdem::accum(dem_dir_fd, mode="d8")
terra::writeRaster(dem_dir_fd, overwrite=T,
"../assets/TIF/dem_chilwa_03_flow_direction.tif")
terra::writeRaster(dem_acc_fd, overwrite=T,
"../assets/TIF/dem_chilwa_04_flow_accumulation.tif")
})
# Compare runtimes
dem_acc_fd = terra::rast("../assets/TIF/dem_chilwa_04_flow_accumulation.tif")
dem_acc_wb = terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation_D8.tif")
cat("\n=== Runtime Comparison ===\n")
NA
NA === Runtime Comparison ===
cat("whitebox flow directions:\n")
NA whitebox flow directions:
cat(sprintf(" Elapsed time: %.3f seconds\n", time_wb["elapsed"]))
NA Elapsed time: 2.023 seconds
cat(sprintf(" User time: %.3f seconds\n", time_wb["user.self"]))
NA User time: 0.016 seconds
cat(sprintf(" System time: %.3f seconds\n\n", time_wb["sys.self"]))
NA System time: 0.013 seconds
cat("flowdem flow directions:\n")
NA flowdem flow directions:
cat(sprintf(" Elapsed time: %.3f seconds\n", time_fd["elapsed"]))
NA Elapsed time: 2.862 seconds
cat(sprintf(" User time: %.3f seconds\n", time_fd["user.self"]))
NA User time: 2.791 seconds
cat(sprintf(" System time: %.3f seconds\n\n", time_fd["sys.self"]))
NA System time: 0.066 seconds
cat("whitebox flow accumulation:\n")
NA whitebox flow accumulation:
cat("Minimum:", min(values(dem_acc_wb), na.rm = TRUE),"cells\n")
NA Minimum: 1 cells
cat("Median:", median(values(dem_acc_wb), na.rm =TRUE), "cells\n")
NA Median: 2 cells
cat("Maximum:", max(values(dem_acc_wb), na.rm= TRUE), "cells\n")
NA Maximum: 1152591 cells
cat("Max contributing area:", round(max(values(dem_acc_wb), na.rm=T) * 0.01, 1), "km²\n")
NA Max contributing area: 11526 km²
cat("flowdem flow accumulation:\n")
NA flowdem flow accumulation:
cat("Minimum:", min(values(dem_acc_fd), na.rm = TRUE),"cells\n")
NA Minimum: 1 cells
cat("Median:", median(values(dem_acc_fd), na.rm =TRUE), "cells\n")
NA Median: 1 cells
cat("Maximum:", max(values(dem_acc_fd), na.rm= TRUE), "cells\n")
NA Maximum: 5017 cells
cat("Max contributing area:", round(max(values(dem_acc_fd), na.rm=T) * 0.01, 1), "km²\n")
NA Max contributing area: 50.2 km²
tmap::tmap_mode("view")
tmap::tm_shape(dem_acc_wb) + tmap::tm_raster(
values="OrRd", title = "Flow Acc.\n(whitebox)",
breaks = c(1, 2, 3, 5, 10, 50, 100, 1000, 1200000),
labels = c("1-2", "2-3", "3-5", "5-10", "10-50", "50-100", "100-1K", ">1K")) +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +
# tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
# tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
# tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_acc_wb
tmap::tm_shape(dem_acc_fd) + tmap::tm_raster(
values="OrRd", title = "Flow Acc. \n(flowdem)",
breaks = c(1, 2, 3, 5, 10, 50, 100, 1000, 1200000),
labels = c("1-2", "2-3", "3-5", "5-10", "10-50", "50-100", "100-1K", ">1K")) +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +
# tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
# tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
# tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_acc_fd
tmap::tmap_arrange(tm_acc_wb, tm_acc_fd, nrow=1)
Figure 2d: Maps comparing flow accumulation results derived from whitebox (left) and flowdem workflows (right).
2.5 Stochastic Depression Analysis
Given the technical challenges of delineating endorheic watersheds without surface outflows, we apply several additional checks below. Specifically, we attempt to characterize depression structures with complementary topographic metrics to account for DEM uncertainty. Specific to this study area, the stochastic approach provides necessary confidence estimates for depression features, which is critical for distinguishing Lake Chilwa’s genuine endorheic basin from spurious artifacts.
Surface depressions in DEMs contain inherent uncertainty from acquisition artifacts, interpolation errors, and measurement precision limitations. The Stochastic Depression Analysis Tool developed in the whitebox library (Lindsay & Creed, 2005) maps topographic depressions while accounting for uncertainty in depression shape resulting from DEM error (Lindsay & Creed, 2006).
The tool uses a Monte Carlo approach where each grid cell follows a Gaussian error probability distribution function (PDF) with mean of zero and standard deviation equal to the LiDAR RMSE. For our 100m DEM derived from multiple sources, we estimate RMSE conservatively at 3m based on vertical accuracy specifications for SRTM data.
We apply 50 iterations following established protocols (Lindsay & Creed, 2005). With each iteration, random samples are drawn from the Gaussian error PDF and added to the original DEM. Depressions in the error-added DEM are filled using Wang and Liu’s efficient algorithm (Wang & Liu, 2006), and affected cells are flagged cumulatively. We classify depressions appearing in ≥80% of iterations as “true depressions” following Lindsay & Creed (Lindsay & Creed, 2006), providing 95% confidence that identified features are not artifacts of DEM error.
# Compute uncertainty using raw and conditioned DEMs
whitebox::wbt_stochastic_depression_analysis(
dem = "dem_chilwa_00_raw.tif",
output = "dem_chilwa_15a_depressions_prob.tif",
wd = "../assets/TIF/",
rmse = 3.0,
range = 10.0,
iterations = 50
)
# Derive binary mask from probability of sinks (~20-80%)
sink_prob = terra::rast("../assets/TIF/dem_chilwa_15a_depressions_prob.tif")
sink_mask_Q1 = sink_prob >= 0.2
sink_mask_Q2 = sink_prob >= 0.4
sink_mask_Q3 = sink_prob >= 0.6
sink_mask_Q4 = sink_prob >= 0.8
sink_mask_Q1[sink_mask_Q1 == 0] <- NA
sink_mask_Q2[sink_mask_Q2 == 0] <- NA
sink_mask_Q3[sink_mask_Q3 == 0] <- NA
sink_mask_Q4[sink_mask_Q4 == 0] <- NA
terra::writeRaster(sink_mask_Q1, overwrite = TRUE,
"../assets/TIF/dem_chilwa_15a_depressions_Q1.tif")
terra::writeRaster(sink_mask_Q2, overwrite = TRUE,
"../assets/TIF/dem_chilwa_15a_depressions_Q2.tif")
terra::writeRaster(sink_mask_Q3, overwrite = TRUE,
"../assets/TIF/dem_chilwa_15a_depressions_Q3.tif")
terra::writeRaster(sink_mask_Q4, overwrite = TRUE,
"../assets/TIF/dem_chilwa_15a_depressions_Q4.tif")
tmap::tmap_mode("plot")
tmap::tm_shape(sink_mask_Q1) + tmap::tm_raster(
values="Greens", title = "Predicted Sinks (>0.2)") +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q1
tmap::tm_shape(sink_mask_Q2) + tmap::tm_raster(
values="Greens", title = "Predicted Sinks (>0.4)") +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q2
tmap::tm_shape(sink_mask_Q3) + tmap::tm_raster(
values="Greens", title = "Predicted Sinks (>0.6)") +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q3
tmap::tm_shape(sink_mask_Q4) + tmap::tm_raster(
values="Greens", title = "Predicted Sinks (>0.8)") +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "turquoise", lwd = 2) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_sink_mask_Q4
tmap::tmap_arrange(
tm_sink_mask_Q4,
tm_sink_mask_Q3,
tm_sink_mask_Q2,
tm_sink_mask_Q1,
nrow=2
)
Figure 2e: Maps comparing thresholds of sink probabilities derived from Monte Carlo simulation of Gaussian distribution (Lindsay & Creed, 2005)
2.5.1 Sink Depth and Wetness
Building on stochastic depression identification, we quantify two complementary metrics:
- Sink depth - vertical distance from each cell to spillover point
- Topographic wetness index (TWI) - water accumulation potential based on upslope contributing area and local slope
Sink depth uses the raw DEM to preserve actual depression depths below spillover points. TWI calculation uses the conditioned DEM and flow accumulation to ensure topologically consistent routing. Together with stochastic mapping, these metrics characterize spatial extent and hydrological characteristics needed to identify terminal drainage points for watershed delineation.
When visually checking results in interactive mapping below, try switching a single layer off using top left drop down button in interactive map windows. Toggling between a stack rasters like this enables quick identification of discrepancies derived from candidate algorithms.
# Calulcate depth of sinks
whitebox::wbt_depth_in_sink(
dem = "dem_chilwa_00_raw.tif",
output = "dem_chilwa_15_sink_depth.tif",
wd = "../assets/TIF/",
zero_background = F
)
# Calculate topographic wetness index
whitebox::wbt_wetness_index(
sca = "dem_chilwa_14_flow_accumulation_D8.tif",
slope = "dem_chilwa_12_filled_wang.tif",
output = "dem_chilwa_16_wetness_index.tif",
wd = "../assets/TIF/"
)
sink_depth <- terra::rast("../assets/TIF/dem_chilwa_15_sink_depth.tif")
wetness_idx <- terra::rast("../assets/TIF/dem_chilwa_16_wetness_index.tif")
flow_accum <- terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation_D8.tif")
tmap::tmap_mode("view")
tmap::tm_shape(sink_depth) + tmap::tm_raster(values = "steelblue",
alpha=1, title = "Sink Depths", style = "quantile", n=7) +
tmap::tm_shape(wetness_idx) + tmap::tm_raster(values="darkgreen",
alpha=1, title="Wetness Index", style = "quantile", n=7) +
tmap::tm_shape(lake) + tmap::tm_borders(col = "turquoise", lwd = 2) +
tmap::tm_layout(legend.text.size = 0.8, legend.title.size = 1)
par(mfrow = c(1, 2))
hist(sink_depth[sink_depth <= 25], breaks = 30,
main = "Depths 0-10m", xlab = "Depth (m)", col = "steelblue")
hist(values(wetness_idx), breaks = 50, xlab = "TWI Value",
main = "Topographic Wetness Index", col = "darkgreen")
Figure 2f: Depression sink depths and their surrounding topographic wetness index scores.
2.6 Terminal Pour Points
Endorheic basins require identification of terminal drainage points representing the ultimate destination of surface flow (Lehner et al., 2022; Zhang et al., 2020) (Li et al., 2025). Unlike exorheic systems with defined outlet pour points, closed basins necessitate locating the convergence point of all inflowing drainage.
We identify terminal points adjacent to maximum flow accumulation within the Lake Chilwa basin. Flow accumulation analysis reveals a critical difference between workflows: whitebox achieves maximum accumulation of 1,152,591 cells (11,526 km²), indicating successful watershed-scale convergence, while flowdem produces only 5,017 cells (50.2 km²), suggesting incomplete drainage network development. This 230-fold discrepancy indicates the flowdem workflow requires troubleshooting before terminal point identification.
Proceeding with the whitebox-derived flow accumulation, we locate terminal points where convergence exceeds 1,000,000 cells, confirming watershed-scale drainage integration (Figure 2g). We test both single-point and multi-point configurations to assess delineation sensitivity to outlet placement in flat terrain. Points are digitized interactively using the flow accumulation raster as reference, positioning outlets at local maxima within the lake basin where multiple drainage pathways converge. The multi-point configuration places four terminal points around the lake perimeter to capture distributed inflow patterns characteristic of endorheic systems with extensive shallow margins.
Terminal points were identified through the following workflow:
- Visual inspection: Examine flow accumulation raster to identify convergence zones exceeding 1,000,000 cells
- Point placement: Use
mapedit::editMap()to enable rapid digitization of pour points at accumulation maxima - Multi-point validation: Add supplementary points around lake perimeter to test sensitivity to outlet configuration
- Reproducibility: Save point locations as shapefiles for workflow replication
The resulting terminal points establish the drainage terminus for watershed delineation, enabling the flowdem::watershed() function to trace all contributing areas converging on Lake Chilwa.
# Inputs
dem_acc = terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation.tif")
streams = dem_acc >= 1000
streams[streams == 0] <- NA
# Manually derive terminal pour points
outlets = mapedit::editMap(mapview::mapView(dem_acc))
outlets = outlets$all |> # convert to sf
sf::st_transform(crs_master) |>
dplyr::select(geometry)
outlets$id <- "terminal_flow"
outlets_add = mapedit::editMap(mapview::mapView(dem_acc))
outlets_add = outlets_add$all |> # convert to sf
sf::st_transform(crs_master) |>
dplyr::select(geometry)
outlets_add$id <- "terminal_flow"
# Outputs
sf::st_write(outlets, "../assets/SHP/outlets.shp", delete_layer=T, quiet=T)
sf::st_write(outlets_add, "../assets/SHP/outlets_multiple.shp", delete_layer=T)
# Visualize
tmap::tmap_mode("view")
tmap::tm_shape(lake) + tmap::tm_borders(col="royalblue") +
tmap::tm_shape(outlets) + tmap::tm_symbols(shape="id",fill="orange", size=0.8) +
tmap::tm_shape(streams) + tmap::tm_raster(values="brewer.reds")NA Reading layer `watershed_chilwa_05_flowdem' from data source `/Users/seamus/repos/map-templates/assets/SHP/watershed_chilwa_05_flowdem.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 1 feature and 1 field
NA Geometry type: POLYGON
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1660000
NA Projected CRS: WGS 84 / Pseudo-Mercator
NA Reading layer `outlets_multiple' from data source
NA `/Users/seamus/repos/map-templates/assets/SHP/outlets_multiple.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 4 features and 1 field
NA Geometry type: POINT
NA Dimension: XY
NA Bounding box: xmin: 3970000 ymin: -1740000 xmax: 3990000 ymax: -1720000
NA Projected CRS: WGS 84 / Pseudo-MercatorFigure 2g: Terminal pour points positioned at maximum flow accumulation within Lake Chilwa basin.
2.7 Watershed Delineation
We delineate the watershed using flowdem::watershed() polygon-based approach, which differs from traditional point-based methods by using the entire lake perimeter as drainage terminus. For spatially extensive terminal sinks such as Lake Chilwa, this function better captures multiple convergent flow paths around the perimeter than a single pour point.
The resulting watershed represents the complete catchment draining to Lake Chilwa through D8-identified surface flow pathways. Post-processing converts the raster watershed output to vector format for integration with other spatial datasets and area calculations. For validation purposes, we present below these flowdem watershed and stream outputs below alongside shapefiles downlaoded externally from OpenTopography collections global datasets.
# Inputs
dem_dir = terra::rast("../assets/TIF/dem_chilwa_03_flow_direction.tif")
dem_acc = terra::rast("../assets/TIF/dem_chilwa_14_flow_accumulation.tif")
streams_test = sf::st_read("../assets/inputs/rivers_site.shp")
watershed_test = sf::st_read("../assets/inputs/watershed_site.shp")
# Delineate flowdem streams
streams = dem_acc >= 1000
streams[streams == 0] <- NA
# Delineate flowdem watershed
watershed = flowdem::watershed(dem_dir, lake) |>
terra::as.polygons() |>
sf::st_as_sf()
streams = terra::crop(
streams, watershed, mask=T)
# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(watershed) + tmap::tm_borders(col = "blue", lwd = 2) +
tmap::tm_shape(streams) + tmap::tm_lines(col = "royalblue", lwd=0.8) +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "royalblue") +
#tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_title("Watershed & Streams (Flowdem)", size=0.8) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_watershed_fd
tmap::tm_shape(watershed_test) + tmap::tm_borders(col = "blue", lwd = 2) +
tmap::tm_shape(streams_test) + tmap::tm_lines(col = "royalblue", lwd=0.8) +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "royalblue") +
#tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
#tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
#tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_title("Watershed & Streams (OpenTopography)", size=0.8) +
tmap::tm_layout(legend.text.size = 0.5, legend.title.size = 0.7) -> tm_watershed_test
tmap::tmap_arrange(tm_watershed_fd, tm_watershed_test, nrow=1)
# Outputs
sf::st_write(watershed, "../assets/SHP/watershed_chilwa_05_flowdem.shp", delete_dsn=T)
terra::writeRaster(streams, "../assets/TIF/streams_chilwa_20a_flowdem.tif", overwrite=T)
#sf::st_write(streams, "../assets/SHP/streams_chilwa_20a_flowdem.shp", delete_dsn=T)
whitebox::wbt_raster_streams_to_vector(
streams = "streams_chilwa_20a_flowdem.tif",
d8_pntr = "dem_chilwa_13_flow_direction_D8.tif",
output = "../assets/SHP/streams_chilwa_20a_flowdem.shp",
wd = "../assets/TIF/")NA Reading layer `streams_chilwa_20a_flowdem' from data source
NA `/Users/seamus/repos/map-templates/assets/SHP/streams_chilwa_20a_flowdem.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 6886 features and 1 field
NA Geometry type: MULTILINESTRING
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1670000
NA Projected CRS: WGS 84 / Pseudo-Mercator
NA Reading layer `rivers_site' from data source
NA `/Users/seamus/repos/map-templates/assets/inputs/rivers_site.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 554 features and 14 fields
NA Geometry type: MULTILINESTRING
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1660000
NA Projected CRS: WGS 84 / Pseudo-Mercator
NA Reading layer `watershed_site' from data source
NA `/Users/seamus/repos/map-templates/assets/inputs/watershed_site.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 1 feature and 16 fields
NA Geometry type: POLYGON
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1660000
NA Projected CRS: WGS 84 / Pseudo-Mercator
NA Reading layer `watershed_chilwa_05_flowdem' from data source `/Users/seamus/repos/map-templates/assets/SHP/watershed_chilwa_05_flowdem.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 1 feature and 1 field
NA Geometry type: POLYGON
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1660000
NA Projected CRS: WGS 84 / Pseudo-Mercator
Figure 2h: Map illustrating watershed delineation derived in flowdem
2.8 Streams Delineation
Stream networks are extracted and classified using hierarchical whitebox processing: (1) extraction by sensitive accumulation threshold, (2) short segment removal, (3) main stem identification, and (4) vectorization.
We apply a sensitive threshold (25 cells) to capture ephemeral channels in marshland, then filter streams shorter than 200m to remove artifacts while preserving primary drainage structure. Main stem identification isolates trunk channels for hydrological modeling applications.
Stream ordering follows Strahler’s method (Strahler, 1957), where first-order streams have no tributaries, and stream order increases when two streams of equal order merge. We retain streams of order ≥3, representing primary drainage channels with contributing areas >10 km².
# Stream Classification
whitebox::wbt_extract_streams(
flow_accum= "dem_chilwa_14_flow_accumulation_D8.tif",
output = "dem_chilwa_18_streams_sensitive.tif",
wd = "../assets/TIF/",
zero_background = T,
threshold = 25,
)
whitebox::wbt_remove_short_streams(
d8_pntr = "dem_chilwa_13_flow_direction_D8.tif",
streams = "dem_chilwa_18_streams_sensitive.tif",
output = "dem_chilwa_19_streams_filtered.tif",
wd = "../assets/TIF/",
min_length = 200
)
whitebox::wbt_strahler_stream_order(
d8_pntr = "dem_chilwa_13_flow_direction_D8.tif",
streams = "dem_chilwa_19_streams_filtered.tif",
output = "dem_chilwa_20_streams_ordered.tif",
wd = "../assets/TIF/",
esri_pntr = FALSE
)
whitebox::wbt_raster_streams_to_vector(
streams = "dem_chilwa_20_streams_ordered.tif",
d8_pntr = "dem_chilwa_13_flow_direction_D8.tif",
output = "../assets/SHP/streams_chilwa.shp",
wd = "../assets/TIF/"
)
# TEST: print(unique(streams_ordered)) >> 8 orders expected
streams_ordered = terra::rast(
"../assets/TIF/dem_chilwa_20_streams_ordered.tif")
streams_major = streams_ordered
streams_major[streams_major < 3] <- 0
streams_major[streams_major == 0] <- NA
terra::writeRaster(streams_major, overwrite = T,
"../assets/TIF/dem_chilwa_21_streams_major.tif")
# Vectorize major streams of interest only
whitebox::wbt_raster_streams_to_vector(
streams = "dem_chilwa_21_streams_major.tif",
d8_pntr = "dem_chilwa_13_flow_direction_D8.tif",
output = "../assets/SHP/dem_chilwa_21_streams_strahler.shp",
wd = "../assets/TIF/"
)
# Tidy
streams_major_sf <- sf::st_read("../assets/SHP/dem_chilwa_21_streams_strahler.shp")
streams_major_sf = st_intersection(streams_major_sf, watershed) |>
dplyr::select("FID") |> dplyr::rename(streams_strahler = FID)
# Outputs
sf::st_write(streams_major_sf, delete_dsn=T,
"../assets/SHP/dem_chilwa_21_streams_strahler.shp")# Inputs
dem = terra::rast("../assets/TIF/dem_chilwa_12_filled_wang.tif")
streams = sf::st_read("../assets/SHP/dem_chilwa_21_streams_strahler.shp")
NA Reading layer `dem_chilwa_21_streams_strahler' from data source `/Users/seamus/repos/map-templates/assets/SHP/dem_chilwa_21_streams_strahler.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 788 features and 1 field
NA Geometry type: MULTILINESTRING
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1670000
NA Projected CRS: WGS 84 / Pseudo-Mercator
watershed = sf::st_read("../assets/SHP/watershed_chilwa_05_flowdem.shp")
NA Reading layer `watershed_chilwa_05_flowdem' from data source `/Users/seamus/repos/map-templates/assets/SHP/watershed_chilwa_05_flowdem.shp'
NA using driver `ESRI Shapefile'
NA Simple feature collection with 1 feature and 1 field
NA Geometry type: POLYGON
NA Dimension: XY
NA Bounding box: xmin: 3910000 ymin: -1790000 xmax: 4040000 ymax: -1660000
NA Projected CRS: WGS 84 / Pseudo-Mercator
wetness = terra::rast("../assets/TIF/dem_chilwa_16_wetness_index.tif")
# Visualize
tmap::tmap_mode("plot")
#tmap::tm_shape(wetness_idx) + tmap::tm_raster(values="darkgreen",
# alpha=1, title="Wetness Index", style = "quantile", n=7) +
tmap::tm_shape(streams) + tmap::tm_lines(col = "blue", lwd=0.8) +
tmap::tm_shape(watershed) + tmap::tm_borders(col = "darkblue", lwd=1.1) +
tmap::tm_shape(lake) + tmap::tm_polygons(fill = "steelblue") +
tmap::tm_layout(legend.text.size = 0.7, legend.title.size = 0.5) +
tmap::tm_scalebar(position=c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark="gray60",position=c("RIGHT", "top")) +
tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_basemap("Esri.WorldTerrain") +
tmap::tm_title("Lake Chilwa's Endorheic Drainage and Stream Network", size=0.8)
Figure 2i: Delineation of Lake Chilwa Drainage Watershed & Classified Stream Network.
2.9 Build 3D map
rivers_strahler = rivers_site |>
dplyr::mutate(
width = as.numeric(
ORD_FLOW
),
width = dplyr::case_when(
width == 3 ~ 16,
width == 4 ~ 14,
width == 5 ~ 12,
width == 6 ~ 10,
width == 7 ~ 6,
TRUE ~ 0
)
) |>
sf::st_as_sf() |>
sf::st_transform(crs = "epsg:4326")
h <- nrow(dem_site)
w <- ncol(dem_site)
dem_matrix = rayshader::raster_to_matrix(dem_site)
dem_matrix |>
rayshader:: height_shade() |>
rayshader::add_overlay(
rayshader::generate_line_overlay(
geometry = rivers_strahler,
extent = dem_site,
heightmap = dem_matrix,
color = "#387B9C",
linewidth = rivers_strahler$width,
data_column_width = "width"
), alphalayer = 1
) |>
rayshader::add_overlay(
rayshader::generate_line_overlay(
geometry = lakes_site,
extent = dem_site,
heightmap = dem_matrix,
color = "#387B9C"
), alphalayer = 1
) |>
rayshader::plot_3d(
dem_matrix,
zscale = 14,
solid = T,
shadow = T,
shadow_darkness = 2,
background = "white",
windowsize = c(600, 600),
zoom = 0.6,
phi = 40,
theta = 0
)2.10 Render 3D map
rayshader::render_highquality(
preview = T,
light = F,
lightdirection = c(135, 45),
lightcolor = c("white"),
lightaltitude = 25,
ambient_light = 0.1,
rotate_env = 0.4,
intensity_env = 0.85,
interactive = F,
parallel = T,
width = w,
height = h,
backgroundhigh="#FFFFFF",
backgroundlow="#FFFFFF"
)
Runtime Log
devtools::session_info()
## ─ Session info ───────────────────────────────────────────
## setting value
## version R version 4.3.0 (2023-04-21)
## os macOS 15.7.1
## system aarch64, darwin20
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Vancouver
## date 2025-11-03
## pandoc 3.8 @ /opt/local/bin/ (via rmarkdown)
## quarto 1.7.33 @ /usr/local/bin/quarto
##
## ─ Packages ───────────────────────────────────────────────
## package * version date (UTC) lib source
## abind * 1.4-8 2024-09-12 [1] CRAN (R 4.3.3)
## backports 1.5.0 2024-05-23 [1] CRAN (R 4.3.3)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.3.3)
## bit 4.6.0 2025-03-06 [1] CRAN (R 4.3.3)
## bit64 4.6.0-1 2025-01-16 [1] CRAN (R 4.3.3)
## bitops 1.0-9 2024-10-03 [1] CRAN (R 4.3.3)
## boot 1.3-32 2025-08-29 [1] CRAN (R 4.3.0)
## brew 1.0-10 2023-12-16 [1] CRAN (R 4.3.3)
## brio 1.1.5 2024-04-24 [1] CRAN (R 4.3.3)
## broom 1.0.8 2025-03-28 [1] CRAN (R 4.3.3)
## bslib * 0.9.0 2025-01-30 [1] CRAN (R 4.3.3)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.3.3)
## callr 3.7.6 2024-03-25 [1] CRAN (R 4.3.3)
## car 3.1-3 2024-09-27 [1] CRAN (R 4.3.3)
## carData 3.0-5 2022-01-06 [1] CRAN (R 4.3.3)
## caret 7.0-1 2024-12-10 [1] CRAN (R 4.3.3)
## class 7.3-23 2025-01-01 [1] CRAN (R 4.3.3)
## classInt 0.4-11 2025-01-08 [1] CRAN (R 4.3.3)
## cli * 3.6.5 2025-04-23 [1] CRAN (R 4.3.3)
## clue 0.3-66 2024-11-13 [1] CRAN (R 4.3.3)
## cluster 2.1.8.1 2025-03-12 [1] CRAN (R 4.3.3)
## codetools 0.2-20 2024-03-31 [1] CRAN (R 4.3.1)
## colorblindcheck * 1.0.2 2023-05-13 [1] CRAN (R 4.3.3)
## colorspace 2.1-1 2024-07-26 [1] CRAN (R 4.3.3)
## cols4all * 0.9 2025-08-28 [1] CRAN (R 4.3.0)
## coro 1.1.0 2024-11-05 [1] CRAN (R 4.3.3)
## countrycode 1.6.1 2025-03-31 [1] CRAN (R 4.3.3)
## covr * 3.6.4 2023-11-09 [1] CRAN (R 4.3.1)
## cowplot * 1.2.0 2025-07-07 [1] CRAN (R 4.3.3)
## crayon 1.5.3 2024-06-20 [1] CRAN (R 4.3.3)
## crosstalk 1.2.2 2025-08-26 [1] CRAN (R 4.3.0)
## curl 7.0.0 2025-08-19 [1] CRAN (R 4.3.0)
## data.table 1.17.8 2025-07-10 [1] CRAN (R 4.3.3)
## DBI 1.2.3 2024-06-02 [1] CRAN (R 4.3.3)
## deldir 2.0-4 2024-02-28 [1] CRAN (R 4.3.3)
## dendextend * 1.19.1 2025-07-15 [1] CRAN (R 4.3.0)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.0)
## DiagrammeR * 1.0.11 2024-02-02 [1] CRAN (R 4.3.1)
## dichromat 2.0-0.1 2022-05-02 [1] CRAN (R 4.3.3)
## digest * 0.6.37 2024-08-19 [1] CRAN (R 4.3.3)
## doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.3.3)
## downlit * 0.4.4 2024-06-10 [1] CRAN (R 4.3.3)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.1)
## dtw * 1.23-1 2022-09-19 [1] CRAN (R 4.3.3)
## dtwclust * 6.0.0 2024-07-23 [1] CRAN (R 4.3.3)
## e1071 1.7-16 2024-09-16 [1] CRAN (R 4.3.3)
## elevatr * 0.99.0 2023-09-12 [1] CRAN (R 4.3.0)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.3)
## evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.3.0)
## exactextractr * 0.10.0 2023-09-20 [1] CRAN (R 4.3.1)
## extrafont 0.19 2023-01-18 [1] CRAN (R 4.3.3)
## extrafontdb 1.0 2012-06-11 [1] CRAN (R 4.3.3)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.3.3)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.3.3)
## flexclust 1.5.0 2025-02-28 [1] CRAN (R 4.3.3)
## flowdem * 0.2 2025-09-14 [1] Github (KennethTM/flowdem@98cdb20)
## FNN * 1.1.4.1 2024-09-22 [1] CRAN (R 4.3.3)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.0)
## foreach 1.5.2 2022-02-02 [1] CRAN (R 4.3.3)
## Formula 1.2-5 2023-02-24 [1] CRAN (R 4.3.3)
## fs 1.6.6 2025-04-12 [1] CRAN (R 4.3.3)
## future * 1.67.0 2025-07-29 [1] CRAN (R 4.3.0)
## future.apply 1.20.0 2025-06-06 [1] CRAN (R 4.3.0)
## gdalUtilities * 1.2.5 2023-08-10 [1] CRAN (R 4.3.0)
## generics 0.1.4 2025-05-09 [1] CRAN (R 4.3.3)
## geodata * 0.6-2 2024-06-10 [1] CRAN (R 4.3.3)
## geojsonsf * 2.0.3 2022-05-30 [1] CRAN (R 4.3.3)
## geos * 0.2.4 2023-11-30 [1] CRAN (R 4.3.3)
## ggmap * 4.0.1 2025-04-07 [1] CRAN (R 4.3.3)
## ggplot2 * 3.5.2 2025-04-09 [1] CRAN (R 4.3.3)
## ggplotify * 0.1.2 2023-08-09 [1] CRAN (R 4.3.0)
## ggpubr * 0.6.1 2025-06-27 [1] CRAN (R 4.3.3)
## ggrepel * 0.9.6 2024-09-07 [1] CRAN (R 4.3.3)
## ggsignif 0.6.4 2022-10-13 [1] CRAN (R 4.3.0)
## ggspatial * 1.1.10 2025-08-24 [1] CRAN (R 4.3.0)
## ggstats * 0.10.0 2025-07-02 [1] CRAN (R 4.3.3)
## giscoR * 0.6.1 2025-08-11 [1] Github (rOpenGov/giscoR@adfed30)
## globals 0.18.0 2025-05-08 [1] CRAN (R 4.3.0)
## glue 1.8.0 2024-09-30 [1] CRAN (R 4.3.3)
## gower 1.0.2 2024-12-17 [1] CRAN (R 4.3.3)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.3)
## gridGraphics 0.5-1 2020-12-13 [1] CRAN (R 4.3.3)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.3.3)
## hardhat 1.4.2 2025-08-20 [1] CRAN (R 4.3.0)
## hdf5r * 1.3.12 2025-01-20 [1] CRAN (R 4.3.3)
## hexbin 1.28.5 2024-11-13 [1] CRAN (R 4.3.3)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.0)
## htmltools * 0.5.8.1 2024-04-04 [1] CRAN (R 4.3.3)
## htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.1)
## httpuv 1.6.16 2025-04-16 [1] CRAN (R 4.3.3)
## httr * 1.4.7 2023-08-15 [1] CRAN (R 4.3.0)
## httr2 * 1.2.1 2025-07-22 [1] CRAN (R 4.3.0)
## interp 1.1-6 2024-01-26 [1] CRAN (R 4.3.3)
## ipred 0.9-15 2024-07-18 [1] CRAN (R 4.3.3)
## iterators 1.0.14 2022-02-05 [1] CRAN (R 4.3.3)
## jpeg 0.1-11 2025-03-21 [1] CRAN (R 4.3.3)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.3.3)
## jsonify 1.2.2 2022-11-09 [1] CRAN (R 4.3.3)
## jsonlite * 2.0.0 2025-03-27 [1] CRAN (R 4.3.3)
## KernSmooth 2.23-26 2025-01-01 [1] CRAN (R 4.3.3)
## knitr 1.50 2025-03-16 [1] CRAN (R 4.3.3)
## later 1.4.4 2025-08-27 [1] CRAN (R 4.3.0)
## lattice * 0.22-7 2025-04-02 [1] CRAN (R 4.3.3)
## latticeExtra 0.6-30 2022-07-04 [1] CRAN (R 4.3.3)
## lava 1.8.1 2025-01-12 [1] CRAN (R 4.3.3)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.3.3)
## leafem * 0.2.5 2025-08-28 [1] CRAN (R 4.3.0)
## leafgl * 0.2.2 2024-11-13 [1] CRAN (R 4.3.3)
## leaflegend 1.2.1 2024-05-09 [1] CRAN (R 4.3.3)
## leaflet * 2.2.2 2024-03-26 [1] CRAN (R 4.3.1)
## leaflet.providers * 2.0.0 2023-10-17 [1] CRAN (R 4.3.3)
## leafpop 0.1.0 2021-05-22 [1] CRAN (R 4.3.0)
## leafsync 0.1.0 2019-03-05 [1] CRAN (R 4.3.0)
## libgeos * 3.11.1-3 2025-03-19 [1] CRAN (R 4.3.3)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.3)
## listenv 0.9.1 2024-01-29 [1] CRAN (R 4.3.3)
## logger 0.4.0 2024-10-22 [1] CRAN (R 4.3.3)
## lubridate * 1.9.4 2024-12-08 [1] CRAN (R 4.3.3)
## luz * 0.5.0 2025-07-29 [1] CRAN (R 4.3.0)
## lwgeom * 0.2-14 2024-02-21 [1] CRAN (R 4.3.1)
## magrittr 2.0.4 2025-09-12 [1] CRAN (R 4.3.0)
## mapedit * 0.7.0 2025-04-20 [1] CRAN (R 4.3.3)
## maptiles * 0.10.0 2025-05-07 [1] CRAN (R 4.3.3)
## mapview * 2.11.2 2023-10-13 [1] CRAN (R 4.3.1)
## MASS 7.3-60.0.1 2024-01-13 [1] CRAN (R 4.3.1)
## Matrix 1.6-5 2024-01-11 [1] CRAN (R 4.3.1)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.3)
## mgcv * 1.9-3 2025-04-04 [1] CRAN (R 4.3.0)
## microbenchmark 1.5.0 2024-09-04 [1] CRAN (R 4.3.3)
## mime 0.13 2025-03-17 [1] CRAN (R 4.3.3)
## miniUI 0.1.2 2025-04-17 [1] CRAN (R 4.3.3)
## ModelMetrics 1.2.2.2 2020-03-17 [1] CRAN (R 4.3.3)
## modeltools 0.2-24 2025-05-02 [1] CRAN (R 4.3.3)
## MPI * 0.1.0 2022-04-05 [1] CRAN (R 4.3.0)
## ncdf4 * 1.24 2025-03-25 [1] CRAN (R 4.3.3)
## nlme * 3.1-168 2025-03-31 [1] CRAN (R 4.3.3)
## nnet * 7.3-20 2025-01-01 [1] CRAN (R 4.3.3)
## openxlsx * 4.2.8 2025-01-25 [1] CRAN (R 4.3.3)
## pacman 0.5.1 2019-03-11 [1] CRAN (R 4.3.3)
## parallelly 1.45.1 2025-07-24 [1] CRAN (R 4.3.0)
## pillar 1.11.0 2025-07-04 [1] CRAN (R 4.3.3)
## pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.3.3)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.3)
## pkgload 1.4.0 2024-06-28 [1] CRAN (R 4.3.3)
## plotly * 4.11.0 2025-06-19 [1] CRAN (R 4.3.3)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.3.3)
## png 0.1-8 2022-11-29 [1] CRAN (R 4.3.3)
## prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.3.3)
## pROC 1.19.0.1 2025-07-31 [1] CRAN (R 4.3.0)
## processx 3.8.6 2025-02-21 [1] CRAN (R 4.3.3)
## prodlim 2025.04.28 2025-04-28 [1] CRAN (R 4.3.3)
## profvis 0.4.0 2024-09-20 [1] CRAN (R 4.3.3)
## progress * 1.2.3 2023-12-06 [1] CRAN (R 4.3.1)
## progressr 0.15.1 2024-11-22 [1] CRAN (R 4.3.3)
## PROJ * 0.6.0 2025-04-03 [1] CRAN (R 4.3.3)
## proj4 * 1.0-15 2025-03-21 [1] CRAN (R 4.3.3)
## promises 1.3.3 2025-05-29 [1] CRAN (R 4.3.3)
## proxy * 0.4-27 2022-06-09 [1] CRAN (R 4.3.3)
## ps 1.9.1 2025-04-12 [1] CRAN (R 4.3.3)
## purrr * 1.1.0 2025-07-10 [1] CRAN (R 4.3.0)
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.3.3)
## randomForest * 4.7-1.2 2024-09-22 [1] CRAN (R 4.3.3)
## rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.3.3)
## raster * 3.6-32 2025-03-28 [1] CRAN (R 4.3.3)
## rasterVis * 0.51.6 2023-11-01 [1] CRAN (R 4.3.3)
## rayshader * 0.37.3 2024-02-21 [1] CRAN (R 4.3.1)
## rayvertex * 0.12.0 2025-02-03 [1] CRAN (R 4.3.3)
## RColorBrewer * 1.1-3 2022-04-03 [1] CRAN (R 4.3.3)
## Rcpp 1.1.0 2025-07-02 [1] CRAN (R 4.3.3)
## RcppParallel 5.1.11-1 2025-08-27 [1] CRAN (R 4.3.0)
## RCurl 1.98-1.17 2025-03-22 [1] CRAN (R 4.3.3)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.3.1)
## recipes 1.3.1 2025-05-21 [1] CRAN (R 4.3.3)
## remotes 2.5.0 2024-03-17 [1] CRAN (R 4.3.3)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.3.0)
## rex 1.2.1 2021-11-26 [1] CRAN (R 4.3.3)
## rgl * 1.3.24 2025-06-25 [1] CRAN (R 4.3.3)
## rgrass * 0.5-2 2025-02-14 [1] CRAN (R 4.3.3)
## rlang 1.1.6 2025-04-11 [1] CRAN (R 4.3.3)
## rmapshaper * 0.5.0 2023-04-11 [1] CRAN (R 4.3.0)
## rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.3.3)
## rpart 4.1.24 2025-01-07 [1] CRAN (R 4.3.3)
## rsconnect * 1.5.1 2025-08-28 [1] CRAN (R 4.3.0)
## RSpectra 0.16-2 2024-07-18 [1] CRAN (R 4.3.3)
## rstatix 0.7.2 2023-02-01 [1] CRAN (R 4.3.0)
## RStoolbox * 1.0.2.1 2025-02-03 [1] CRAN (R 4.3.3)
## rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.3.3)
## rts * 1.1-14 2023-10-01 [1] CRAN (R 4.3.3)
## Rttf2pt1 1.3.12 2023-01-22 [1] CRAN (R 4.3.3)
## s2 * 1.1.9 2025-05-23 [1] CRAN (R 4.3.3)
## sass 0.4.10 2025-04-11 [1] CRAN (R 4.3.3)
## satellite 1.0.6 2025-08-21 [1] CRAN (R 4.3.0)
## scales * 1.4.0 2025-04-24 [1] CRAN (R 4.3.3)
## sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.3.3)
## sf * 1.0-22 2025-08-25 [1] Github (r-spatial/sf@3660edf)
## shiny 1.11.1 2025-07-03 [1] CRAN (R 4.3.3)
## shinyjs 2.1.0 2021-12-23 [1] CRAN (R 4.3.0)
## shinyWidgets 0.9.0 2025-02-21 [1] CRAN (R 4.3.3)
## slippymath 0.3.1 2019-06-28 [1] CRAN (R 4.3.0)
## sp * 2.2-0 2025-02-01 [1] CRAN (R 4.3.3)
## spacesXYZ 1.6-0 2025-06-06 [1] CRAN (R 4.3.3)
## spData * 2.3.4 2025-01-08 [1] CRAN (R 4.3.3)
## spdep * 1.4-1 2025-08-31 [1] CRAN (R 4.3.0)
## stars * 0.6-8 2025-02-01 [1] CRAN (R 4.3.3)
## stringi 1.8.7 2025-03-27 [1] CRAN (R 4.3.3)
## stringr * 1.5.2 2025-09-08 [1] CRAN (R 4.3.0)
## supercells * 1.0.0 2024-02-11 [1] CRAN (R 4.3.1)
## survival 3.8-3 2024-12-17 [1] CRAN (R 4.3.3)
## svglite 2.2.1 2025-05-12 [1] CRAN (R 4.3.3)
## systemfonts 1.2.3 2025-04-30 [1] CRAN (R 4.3.3)
## terra * 1.8-60 2025-07-21 [1] CRAN (R 4.3.0)
## terrainr * 0.7.6 2025-07-25 [1] CRAN (R 4.3.0)
## testthat * 3.2.3 2025-01-13 [1] CRAN (R 4.3.3)
## textshaping 1.0.3 2025-09-02 [1] CRAN (R 4.3.0)
## tibble * 3.3.0 2025-06-08 [1] CRAN (R 4.3.3)
## tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.3.1)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.3.1)
## tidyterra * 0.7.2 2025-04-14 [1] CRAN (R 4.3.3)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.0)
## timechange 0.3.0 2024-01-18 [1] CRAN (R 4.3.3)
## timeDate 4041.110 2024-09-22 [1] CRAN (R 4.3.3)
## tmap * 4.1 2025-05-12 [1] CRAN (R 4.3.3)
## tmaptools * 3.2 2025-01-13 [1] CRAN (R 4.3.3)
## torch 0.16.0 2025-08-21 [1] CRAN (R 4.3.0)
## traudem * 1.0.3 2025-10-28 [1] Github (lucarraro/traudem@a30c9c9)
## tzdb 0.5.0 2025-03-15 [1] CRAN (R 4.3.3)
## unifir 0.2.4 2024-02-01 [1] CRAN (R 4.3.3)
## units 0.8-7 2025-03-11 [1] CRAN (R 4.3.3)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.3)
## usethis 3.1.0 2024-11-26 [1] CRAN (R 4.3.3)
## uuid 1.2-1 2024-07-29 [1] CRAN (R 4.3.3)
## V8 6.0.4 2025-06-04 [1] CRAN (R 4.3.3)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.3)
## viridis 0.6.5 2024-01-29 [1] CRAN (R 4.3.1)
## viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.3)
## visNetwork 2.1.2 2022-09-29 [1] CRAN (R 4.3.0)
## whitebox * 2.4.3 2025-10-28 [1] Github (opengeos/whiteboxR@294ea4b)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.3.3)
## wk 0.9.4 2024-10-11 [1] CRAN (R 4.3.3)
## xfun 0.53 2025-08-19 [1] CRAN (R 4.3.0)
## xgboost * 1.7.11.1 2025-05-15 [1] CRAN (R 4.3.3)
## XML 3.99-0.18 2025-01-01 [1] CRAN (R 4.3.3)
## xml2 1.4.0 2025-08-20 [1] CRAN (R 4.3.0)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.3)
## xts * 0.14.1 2024-10-15 [1] CRAN (R 4.3.3)
## yaml 2.3.10 2024-07-26 [1] CRAN (R 4.3.3)
## yulab.utils 0.2.1 2025-08-19 [1] CRAN (R 4.3.0)
## zeallot 0.2.0 2025-05-27 [1] CRAN (R 4.3.3)
## zip 2.3.3 2025-05-13 [1] CRAN (R 4.3.3)
## zoo * 1.8-14 2025-04-10 [1] CRAN (R 4.3.3)
##
## [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
## * ── Packages attached to the search path.
##
## ──────────────────────────────────────────────────────────